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Abstract 

Recently, the DSBGK method {note: the original name DS-BGK is changed 
to DSBGK for simplicity) was proposed based on the BGK equation to re- 
duce the stochastic noise in simulating rarefied gas fiows at low velocity, 
in which the deviation from equilibrium state is small making the tradi- 
tional DSMC simulation time-consuming due to the dominance of noise in 
transient results. In both DSMC and DSBGK simulations, the simulated 
molecules move into and out of cells randomly and frequently. Consequently, 
the transient information of molecules in each particular cell contains sig- 
nificant noise. The DSMC method uses the transient values of molecular 
variables to compute the cell's variables (including number density, fiow ve- 
locity and temperature) and so the stochastic noise in its cell's variables is 
remarkable particularly in the case of low velocity. In the DSBGK simulation, 
the increments rather than the transient information of molecular variables 
are used to update the cell's variables based on the mass, momentum and 
energy conservation principles of intermolecular collision process. This up- 
dating scheme significantly reduces the noise in cell's variables of DSBGK 
simulations because the molecular variables are updated smoothly by the ex- 
trapolation of acceptance-rejection scheme and so their increments contain 
low noise. The detailed comparisons of algorithms and results between the 
DSMC and DSBGK methods are given here. Several benchmark problems 
are simulated to verify the DSBGK method by comparison with the DSMC 
method as criterion. 
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1. Introduction 

For micro gas fiows, the Boltzmann equation rather than the Navier- 
Stokes equation should be used due to high Knudsen number Kn = X/L 
where A is the molecular mean free path and L is the characteristic length 
of the fiow problem. In addition, the infiuence of boundary condition to the 
solutions becomes dominant because the frequency of molecular refiection on 
the solid wall, compared to the frequency of intermolecular collision, increases 
with Kn. Unfortunately, the characteristic velocity of micro gas fiows is usu- 
ally much smaller than the molecular random thermal velocity and sometimes 
the variations of quantities of interest inside the fiow domain are very small, 
which makes the traditional DSMC method [1] time-consuming although it 
is successful in the case of high velocity. 

The DSBGK method [2j was proposed to improve the eflBciency in sim- 
ulating micro gas fiows and verified in the lid-driven, Couette and channel 
fiow problems [^-0 by comparison with the DSMC method as criteria. The- 
oretically, it can be proved, as will be discussed later, that the solution of the 
DSBGK method converges to the steady-state solution of the BGK equation 
[4j. The application of the CLL refiection model [5j-[6j in the DSBGK method 
is possible and few tentative results were compared with the DSMC results in 
[5. Although based on the BGK equation obtained by using a simple model 
to replace the intermolecular collision integral of the Boltzmann equation, 
the DSBGK method agrees well with the DSMC method at Kn = 0.063 and 
6.3 in the lid-driven problem [2]. This is because the molecular refiection 
on wall, the dominant eflFect in micro gas fiows, is modeled by the DSBGK 
method in the same way as by the DSMC method. Theoretically, the er- 
ror due to simplification to the intermolecular collision process vanishes and 
the solution depends only on the boundary condition when Kn 00. The 
DSBGK method achieves high efiiciency by avoiding generating random frac- 
tions in the intermolecular collision process and using the increments (instead 
of transient values) of molecular variables to update cell's macro quantities, 
which significantly reduces the statistical noise due to discontinuous events 
of simulated molecules moving into and out of cells. Consequently, the to- 
tal computational time used by the DSBGK simulation almost not increase 
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with the decrease of magnitude of the deviation from equihbrium state and 
sometimes the average process can be avoided as the transient ceU's variables 
contain few stochastic errors [2j. In addition to its high-efficiency, the DS- 
BGK method has many numerical advantages including simplicity, stability, 
convenience for complex configuration and for parallel computation because 
the basic algorithmic structure of the DSMC method is employed. 

The comparison between the DSMC and DSBGK algorithms is given 
here. Theoretical analysis is provided to show the convergence of the DS- 
BGK method to the BGK equation. Then, the results of several benchmark 
problems, including the Couette fiow, channel fiow, lid-driven fiow and ther- 
mal transpiration problem, are listed together to show the agreement of the 
DSBGK method with the DSMC method. The benchmark problems are di- 
vided into closed and open problems to discuss the efficiency and stability 
of the DSBGK simulation separately. In closed problems, the long-period 
ffuctuation is observed in the number density distribution of DSBGK simu- 
lations. Many simulated molecules are employed to reduce the magnitude of 
ffuctuation and improve the numerical stability. Consequently, the memory 
usage is increased remarkably but the efficiency is still very high as shown 
in the closed lid-driven problem [2j. In open problems, the boundary con- 
dition with ffxed number density eliminates the unphysical ffuctuation and 
the DSBGK simulation remains stable even when using about 10 simulated 
molecules per cell, which signiffcantly reduces the memory usage and so im- 
proves the applicability in open problems of large scale. 

2. DSMC Method 

The DSMC method [Ij, which is successful in simulating rareffed gas 
ffows at high velocity, was proposed based on physical understanding with 
appropriate theoretical analysis. In fact, the DSMC algorithm in simple 
cases can be understood by using the importance sampling scheme to solve 
the Boltzmann equation ^-^Sj, which is discussed here. The rareffed gas 
ffow is described by the Boltzmann equation. We consider gas ffows of single 
component in the absence of external body force. If the molecule is modeled 
by a hard sphere with ffxed diameter D, the Boltzmann equation is: 



dt ^ dxj dt 
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/(t, X, c) is the unknown probability distribution function, t is the time, x is 
the spatial coordinate and c is the molecular velocity, /i = /(t, x, Ci) and /2 = 
/(t, X, C2), the delta function 5i 5(ci — c), (^2 ^^(^2 — c), (^(c^^ — c), ^2 = 
5(0^2 — c), 5^ = |c2 — Ci|, the post-collision velocities c^i,c^2 determined by 
the pre-collision velocities Ci, C2 and the solid angle f], df] = ^iiiLpdLpdO where 
(/::' G [0, tt] is the polar angle (the deflection angle in intermolecular collisions) 
and 6 G [0, 27r] is the azimuthal angle of the spherical coordinate system. The 
total collision section is = J^^ D'^/AdQ = ttD'^. The boundary condition 
will be discussed later in section [3^ together with the DSBGK method. After 
getting the solution of c), the number density n(t,x), flow velocity 

u{t^ x) and temperature T(t, x) are computed 
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where m is the molecular mass and /cb is the Boltzmann constant. Higher 
order momentums, like shear stress tensor and heat flux, are computed sim- 
ilarly. 

In the DSMC simulation [Ij, each simulated molecule / carries two molec- 
ular variables: position xi and velocity q. In order to reduce the memory 
usage, the number of simulated molecules is much smaller than that of the 
real molecules contained in the flow domain and so we assume that each 
simulated molecule represents number of real molecules. Note that N 
is a constant and very large to make each cell usually containing about 20 
simulated molecules. The molecular position and velocity are selected at the 
initial state and updated during the simulation process appropriately such 
that the set of all simulated molecules [x/,Q]aii represents the probability 
distribution function / and its evolution with time, which means that the 
simulated molecules are distributed according to / in the phase space (x, c) 
at any moment t. The flow domain is divided into many cells and n^u^T are 
estimated by summation inside each cell k using N /Vk to replace /dc in Eq. 
([2]) as fdcdx is the number of real molecules in the velocity space element 
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dc and the physical space element dx: 
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where Vk is the volume of ceU fc, ^ is the summation over those simulated 
molecules located inside cell k at any particular moment t. For example, 
is the product of and the number of simulated molecules and so 
equal to the number of real molecules inside cell k. 

df . ^ a/, df 



During each time step At, we split — into — I 

at at 

df 



due to free 



molecular motions and — |coii due to intermolecular collisions. As [x/,Q]aii 



dt 

is a representative sample of /, 
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3 is represented by updating xi when 

simulated molecules move uniformly and in a straight line. 
df 

For — Icoib we need to calculate the increment A/Loii of / after each At 
dt ^ 

at all spatial points x and all velocity points c inside the whole phase space. 
In order to make A/|coii tractable, we assume that the coordinates xi of those 
simulated molecules inside the same cell k are the same (notated by x center, k)- 
Then, we only need to compute A/|coii at those discrete spatial points Xcenter.k 
of each cell (as / = and so A/ = at other spatial points) but still at all 
velocity points. The distribution function at Xcenter,/c is fk = ^^{ci — c)N/Vk 
which is consistent with Eqs. ([2])-([3]) as Uk = /^dc = ^A^/V/e (again, 
^ is over simulated molecules inside cell k). At the end of each At and for 
each cell fc, we compute Afk\co\\ according to the Boltzmann equation: 
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where M = -— and G = —(5^ + S[ -82- 5i)- r , fk^i = 

^ 5(q — Ci)N/Vk is the distribution function of Ci at Xcenter,/c- Note that the 

value of g has upper bound here as fk is nonzero only at finitely many discrete 

velocity points q. Although the value of (5'crT)max can be any constant in Eq. 

([4]) , it should be updated appropriately by the existing values gaT in all cells 

during each At in the DSMC simulation such that the ratio g(TT/{g(^T)ina^ 

is always {note: practically will be 'almost always') smaller than 1 which is 

required by the following acceptance-rejection scheme. But, if (5'crT)max is 

much larger than that required to make all ratios smaller than 1, the number 

M of tentative collision pairs is very large making the simulation process 

time-consuming due to low acceptance probabilities of the tentative collisions 

4 fk 1 

(see the following analysis). Note that - — dQ = 1, — -dci = 1, 
^ & J ; Jo 4^^ ' J-00 

fk 2 

— -dc2 = 1 and so, Afk\co\\ is equal to M < G > where < G > is the 
^ rik 

expected value of G. We use the importance sampling scheme to estimate 
< G >, namely EsampieG^j ^< ^ > where Esampie is the sum 

' ^sample 

of ^sample Humbcr of representative Gj = Gj(f],Ci,C2) with f],Ci,C2 being 

D f]^ f]^ 2 

selected according to their probability densities - — , — — respectively. 

4crT rik rik 

Furthermore, we let risampie = M and so Afk\co\\ ^ Esampie 

For any Gj, we select particle ji randomly and uniformly from those 
simulated molecules inside cell k and thus Ci = Cj^ is selected according to 
fk.i/'^k as required because fk^i = E^(^^ ~ Ci)A^/V/e, which implies that all 
simulated molecules should be selected equivalently. The number of simu- 
lated molecules inside dc represents fk- Then, we select particle ^2 (^2 7^ ji) 
randomly and uniformly inside cell k and use Cj^ as the representative 
value of C2, which also implies that C2 is selected according to fk,2/'^k where 
fk,2 = E^(^^ ~ C2)N/Vk- As i^^/(4crT) is a constant, we select Q randomly 
and uniformly from the whole surface of unit sphere, which is equivalent 
to selecting the post-collision , c^-^ randomly by the hard-sphere collision 
model as Q is used only to calculate , . Now, we have Cj^ , Cj^ , , and 
9j = ~ ^j2\- Assuming that (5'crT)max is known for the current At, Gj is 

equal to 77 [^^(^2 - + S{^j^ - c) - S{cj^ - c) - S{cj^ - c)]—^^^ — . Now, 
the acceptance-rejection scheme is used to handle the fraction '^•^^^ — . If 

(^CrT)max 
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^■^ ^ > Rf where Rf is a random fraction distributed uniformly inside 



N 

[0, 1] , we let Gj = — [5(4 -c) + S{^j^ -c)- S{cj, - c) - S{cj, - c)] and Gj = 

otherwise. Note that fk = I]5(q - c)N/Vk and Afk\co\\ - Zlsampie^j 

so fk becomes ^(^(q — c)N/Vk + Ssampie^j ^fter intermolecular coUisions. 

This imphes that if - — ^— > i?/, — [(^(cj2 — c) + (5(cj^ — c)] contained in 

^ 7^^(q — c) is canceled by — [— (^(cj2 — c) — S{cj^ — c)] contained in Gj and 
Vk Vk 

N N ^ 

meanwhile Ty[S{c^j^—c) + 5{^- —c)] contained in Gj is added to ^ — (^(q — c), 
Vk Vk 
N N 
namely replacing —[5{cj^ - c) + 5{cj-, - c)] by TrlK^n ~ ^ + ^(^ji ^] 

»^/c ^k 

N 

^ — 5(q — c). Till now, the replacement may contribute nothing if we are 
Vk 

discussing A//e|coii at velocity points c different from Cj^, Cj2, because 

N ^ _ N 

both the original — [^(^73 — c)+5(co^ — c)] and the new — [5(c^- —c)-\-5{^- — c)] 

are equal to zero at those c. So, we consider A//e|coii at all velocity points 
c together and specify that the same set of samples Gj is used to compute 

A//e|coii at all different c. Then, if ^"^^^ — > i?/, the contribution of Gj 

to A//e|coii in the whole velocity space is nonzero only at four velocity points 
and equivalent to changing the velocities Cj^ , Cj^ to c^-^ , , respectively, which 
means that a pairwise intermolecular collision happens. So, we select M num- 
ber of tentative collision pairs for each cell k at the end of each At and use 
^^^^ of each pair ^1,^2 as the acceptance probability to judge whether a 



(^^T)r 

pairwise collision happens. This is the algorithm used in the DSMC method. 

For dense fluids, the importance sampling scheme was used in [9j to solve 
the Enskog equation, which is an extension of the Boltzmann equation by 
considering the intermolecular repulsive force at short distance but still ne- 
glecting the intermolecular cohesive force at long distance. The cohesive force 
is vital in simulating two-phase flows [TO]. For problems at low velocity, the 
intermolecular collision integral of the Boltzmann equation is simplified and 
evaluated by the importance sampling scheme to improve the efficiency in 
the LVDSMC method [11], which conserve the mass on average. A scheme 
was proposed in [12j to conserve the mass strictly. 
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3. DSBGK Method 



We consider gas flows of single component. In the absence of external 
body force, the BGK equation [4J can be written as a Lagrangian form: 



d/^a/ df 
dt dt ^^dx 



where /(t, x, c) is the unknown probability distribution function, t is the time, 
X is the spatial coordinate and c is the molecular velocity, the parameter 
V is selected appropriately to satisfy the coeflicient of viscosity /x or heat 
conduction [13j : 

Mbgk — 

5kB nk^T 

^BGK = 

Im V 

and the Maxwell distribution function /gq is: 

/eq t,x,c) = n( — -—f/^exp \ rr. ] 7 

where m is the molecular mass and /cb is the Boltzmann constant, the number 
density n, flow velocity u and temperature T are functions of t and x and 
deflned by Eq. ^ using / . 

In the DSBGK method [2j, the simulation process is divided into a se- 
ries of time steps At and the flow domain is divided into many cells. The 
selections of At and cell size are the same as in the DSMC method when 
simulating problems of high Kn. Many simulated molecules are employed 
to represent the distribution function / and its evolution with time. The 
main idea of this method is to track down the evolution of / along enormous 
molecular trajectories at constant velocities, which are selected randomly 
when simulated molecules are generated or reflected at the boundaries. Each 
simulated molecules / carries four molecular variables: position x/, molec- 
ular velocity q, number Ni of real molecules represented by the simulated 
molecule /, and Fi which is equal to the representative value f{t^xi^ci) of / 
at the moment t and point (x/,q) in the phase space, [x^, q, A^/]aii is a (not 
unique) representative sample of / and [F/]aii is the representative value of /. 
The compatibility condition, namely [x^, q, A^/]aii and [F/]aii are related to the 
same /, is required during the simulation process. Note that the evolution 
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of / is due to three factors: free molecular motion, intermolecular collision 
and molecular reflection on the wall. 

For the evolution of / due to free molecular motions and intermolecular 
collisions, [F/]aii is changed and then [x/, q, A^/]aii is updated correspondingly 
by changing xi and Ni rather than q. Note that Ni is a constant and xi is 
changed alone to represent the evolution of / due to free molecular motions 
and then q is changed randomly to represent the evolution of / due to in- 
termolecular collisions in the DSMC simulation. The position xi is updated 
along the trajectory of molecular free motion. Fi is updated with xi by the 
Lagrangian description of the BGK equation where /eq is replaced inside each 

ceU k by the transitional /eq,tr,/c = ^tr,fe( ^ , ^ )^/^exp[ ^ — ]. 

The cell's variables ntr,)^^ '^tr,^, ^tr,/c are updated by x/,q and the increment 
(instead of transient value) of Ni based on the mass, momentum and energy 
conservation principles of intermolecular collision process. Note that we use 
the subscript tr to distinguish the transitional cell's variables ntr,^, ^t^.k^ T^r^k 
from rik.Uk.Tk, which are computed by the transient x^, q, Ni as in Eq. ^ 
because [x/, q, A/'/Jaii is a representative sample of /. The increment of Ni is 
due to the intermolecular collision effect and computed by the extrapolation 
of acceptance-rejection scheme, which avoids the time-consuming process of 
frequently generating random fractions and employs the changing informa- 
tion of Fi making the compatibility condition satisfled. 

For the evolution of / due to molecular reflection at xi on the wall, q 
is changed to Q,new = Cr + ixwaii where is the reflecting velocity selected 
randomly in the local Cartesian reference system moving at the wall velocity 
iTwaii- But, A^^ remains unchanged to conserve mass. Then, Fi is updated to 
^/,new = fit-i^u^i^new)-) which also satisflcs the compatibility condition. 

3.1. Initialization process 

At the initial state, the cell variables ntr,/c, '^tr,fc, T^^.k are equal to the initial 
macro quantities which are usually uniform. The initial molecular position 
xi and velocity q are selected randomly as in the DSMC simulation and then 
Fi is equal to /eq,tr,iEc(0, xi^ci). The initial values A^^^^^o of A'^ for different sim- 
ulated molecules are usually the same and selected appropriately such that 
the total number of simulated molecules, which is equal to A^totai,reai/^/,t=o 
where A^totai,reai is the total number of real molecules, takes a acceptable value. 
The smaller the value of A^^^^^o is, the larger the total number of simulated 
molecules at the initial state will be. 
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3.2. Algorithms for molecular motion and intermolecular collision 



A particular molecular trajectory during one time step 
is divided into three segments by the cell's interfaces 




Summation over segments (blue color) in the same cell 
(different lengths of trajectories due to different velocities) 




(a) Trajectory division 



(b) Summation inside each cell 



Figure 1: Schematic models of the DSBGK simulation. 



In the DSBGK simulation, each simulated molecule moves uniformly and 
in a straight line before encountering the boundary. As we can see from 
Fig. [T} the molecular trajectory during each At may be divided into several 
segments by cell's interfaces or remains as a single segment if not yet arriving 
at any cell's interface at the end of current At. As each segment is located 
inside a particular cell fc, is conveniently updated along each segment in 
sequence according to the Lagrangian form of the BGK equation using /eq,tr,fc 
of cell k. Note that /eq,tr,/c is constant for a particular simulated molecule / 
and cell k as ntr,^:, '^tr,^, ^tr,fc and q are fixed. So, Fi is updated using Eq. ([s]) 

dFi 

obtained by finishing the integration of Eq. ([5]), namely — — = '^(/eq,tr,fc — ^/), 
with respect to t along the segment concerned: 

Fl^new = /eq,tr,fc + {Fl - /eq,tr,/c) eXl[){-vAkti) (8) 

where Fi is the previous value and F^^new is the new value after the inter- 
molecular collision, A^ti is the time interval used by the simulated molecule 
/ during the current At to go through the segment inside cell k. As the 
molecular trajectory is divided first by the time step At and then by the 
cell's interfaces, A^ti < At. If the trajectory during the current At is di- 
vided into several segments by the cell's interfaces, A^ti < At and Eq. ([s]) is 
used repeatedly to update Fi for the consecutive segments in sequence. After 
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updating Fi for each segment, Ni is updated correspondingly: 

A^/,new = NiFi^^,^/Fi (9) 

which is based on the extrapolation [2j of acceptance-rejection scheme that if 
[xu Q, A^^]aii is a representative sample of /i, [xu Q, A^K/2//i)/]aii is a represen- 
tative sample of /2, where {f2/ fi)i is the ratio of /2 and /i at the same point 
(t^xi^ci). Equation ^ could be understood by considering two steps: in 
the first step without intermolecular collision, xi is updated with t along the 
trajectories but q, F/, Ni keep unchanged as f{t + At, x + Ate, c)=f{t^ x, c); 
then, F/ is changed to Fi^new due to intermolecular collision and t^xi^ci keep 
unchanged, so, Ni is changed correspondingly to A^/,new by Eq. The 
precondition of using the extrapolation of acceptance-rejection scheme is 
that [x/, Q, A^/]aii is a representative sample of / whose representative value 
is [Fijaii before intermolecular collision, namely the compatibility condition 
must holds before using the extrapolation of acceptance-rejection scheme. 
Then, the updating algorithms of x/, F/, Ni with t for the free molecular mo- 
tion and intermolecular collision processes make the compatibility condition 
constantly satisfied due to using the extrapolation of acceptance-rejection 
scheme. In the molecular refiection process on the wall, the compatibility 
condition is satisfied automatically. 

The idea of the updating algorithms along molecular trajectories at con- 
stant velocities is inspired by the Lattice Boltzmann method (LBM). In turn, 
the physical understanding of the kinetic equation is also helpful to the de- 
velopment of LBM algorithm. Recently, an alternative scheme was proposed 
in [llj to compute the strain rate tensor for the application of large eddy 
simulation (LES) in the LBM. 

The cell's variables Utr^k^ '^tr,/c5 rtr,/c are used in Eq. ^ to determine /eq,tr,fc 
and updated at the end of each At. During the current At and for each 
cell k (see Fig. [l] right), some simulated molecules run inside cell k and 
their increments Aj^Ni = A^/,new — inside cell k are already known. A^Ni 
is the number increment of real molecules of class q associated with the 
intermolecular collisions inside cell k during the current time step. We make 
summation ^ AkNi over those simulated molecules running inside cell k 
during the current At {note: simulated molecule / may contribute more 
than one term to the summation if it refiects on the wall back into the cell 
k). Note that A^^A^^ in this summation is the increment information rather 
than transient information in the summation of Eq. ^ used in the DSMC 
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method. Obviously, ^ A^Ni means the number increment of real molecules 
of all existing classes associated with the intermolecular collisions inside the 
same cell k during the same time step. So, ^ ^k^i is expected to be zero as 
required by the mass conservation principle. Usually, this summation is not 
exactly equal to zero due to numerical error. So, we decrease ntr,/c if ^ ^k^i 
is positive and then ^ A^Ni will decrease at the next At as each term A^Ni 
decreases due to Eqs. ([8])-([9]), and vice versa. It works as an auto-regulation 
scheme which makes "^Aj^Ni approaching to zero. Similarly, "^(AkNimci) 
and ^{AkNimc^ /2) are related respectively to the momentum increment and 
kinetic energy increment of real molecules of all existing classes associated 
with the intermolecular collisions inside the same cell k and during the same 
At. They are expected to be zero according to the momentum and energy 
conservation principles of intermolecular collision process and so can be used 
to update 'Utr,^ and T^r^k by auto-regulation schemes. The auto-regulation 
schemes are: 

' new ntr^kVk-Y.^kNl 



tr,/c T/ 



U 



k 

ntr,kVkUtr,k - Y^i^kNiCi) 



'hr^k ^k 



Tl 



new 



[ntv,kVk[ 7, \ n — > ~ l^K^k^yi^^)] - n^r,k ^ 



'hr,k^k 2 

(10) 

where n^^'^^u^^'^^T^^^ are the new values of number density ntr,^, flow ve- 
locity 'Utr,/c and temperature T^j-^k of cell fc, respectively, Vk is the volume 
of cell k. The updating schemes of Eq. (10) make J^^k^^h ^(^a^^/q), 
^{AkNim^ /2) converging to zero and then nt^^k->^tr,k->Ttr,k will fluctuate 
around their steady state solutions due to stochastic eflFect. 

We use A^^, Fi to represent the previous values at the origin of the segment 
located inside cell k during the current At and use Ni^^^^^ Fi^^^^ for the new 
values at the end of that segment after intermolecular collision as in Eqs. 
([8])-([9]). Note that any possible representative trajectory is selected accord- 
ing to its probability (see section 3.4) as the molecular reflecting velocity is 
selected randomly according to the boundary reflection model. Thus, it can 
be expected that the feature of all existing classes represents the feature of all 
possible classes and so the summation over all existing classes is equivalent 
to the integration over all possible classes like replacing Eq. ^ by Eq. ^ 
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in the DSMC simulation. We replace Ni by 14F/dQ where dc/ is the velocity 
space element around q as the compatibility condition is satisfied. Note that 
Akti is the time interval used by the simulated molecule / inside cell k dur- 
ing the current At and so A^ti = At for those simulated molecules moving 
inside the same cell (namely the trajectory during the current At is a single 
segment without division by the cell's interfaces). We assume that At is very 
small making most simulated molecules moving inside the same cell during 
each At and so A^ti At. The integral expression of mass conservation of 
the DSBGK simulation for each cell k is: 

= Y^iVkFi^ne^dCi) - Y^iVkFldCi) 
= ^[^^(/eq,tr,fc " Fi)AktldCl] 

/OO 
v{feq,tT,k - /)dc 
-OO 

dFi 

where = '^(/eq,tr,A; " Fi) and the last approximate equality holds as 
Akti At and Fi is the representative value of /. So, after convergence 
with ^ AkNi — 0, we get 

/»oo 

^(/eq,tr,fc-/)dc = (12) 



/ 

J — c 



But, it is not necessary to require At being very small. If At is large, 
^ Aj^Ni = still implies ^(/eq,tr,/c — /)dc = as both of them represent 
the mass conservation of intermolecular collision process of the same evolu- 

tion equation — = v{feq,tr,k — /)• After convergence with ^(A/^A^/q) = 

and ^(A/jA^/mcf /2) = 0, the integral expressions of momentum and energy 
conservations can be obtained similarly. So, the following equation is satisfied 
for each cell k after convergence: 



/oo 
-oo 



^'(/eq,tr,fc - /)^idc = (13) 

I 

where -01 = 1, {i'2,i'3,'ip4) = c^ij)^ = m(?/2. As the original BGK equation 
satisfies j^^v{f^ - f)tpidc^ 0, we have X!^(/eq,tr,fc - /eq)^idc = which 
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implies that ritr^k = ^5 '^tr,/c = '^5 ^tr,/c = T for each ceU after convergence 
according to the definitions of /eq,tr,/c and /eq. 

So, the solutions of ntr,/c, '??tr,fc, Ttr,^ of the DSBGK method are the dis- 
crete solutions of n, T of the BGK equation after convergence under the 
same boundary condition. Then, the transitional feq,tT,k used in the DSBGK 
method is equal to the original /eq of the BGK equation inside each cell k. 
Consequently, [F^jaii and [x^, q, A^^jaii are the representative value and sample, 
respectively, of the solution / of the BGK equation, which implies that any 
higher-order moment, including stress tensor and heat fiux, calculated by 
the DSBGK method agrees with that obtained by solving the BGK equation 
using other numerical methods as in [15] among others. 



Note that the updating scheme of Eq. (10) conserve the 'total' value 
ntr,kVk + Y.^i inside each ceU k as {n^^J^Vk - ritr^kVk) + ^k^i = {note: 
'total' with quotation marks means the sum of cell quantity and molecular 
quantity). So, ^Domain ^tT,kVk + Y.Dom^in coustaut during the simulation 
process because ntr,/c and Ni are unchanged during the molecular refiection 
process on the wall {note: the summation ^oo^iain ^^^^ whole fiow 
domain, namely over all cells and all simulated molecules, respectively). The 
'total' momentum and energy of simulated molecules and cells are unchanged 
when using Eq. (10) but not conserved during the whole simulation process 
due to molecular refiections on the wall, which conserve the mass but not 
momentum and energy. Note that the conservations of the 'total' mass, 
momentum and energy by the updating scheme of Eq. (10) are artificial 



restrictions. Eq. ([10|) can be modified by adding arbitrary diflFerent positive 
factors before ^ A/^A^/, ^(A/^A^/q), ^(Aj^A^/mcf /2) to regulate the conver- 
gence speed in open problems. But, the 'total' mass should be conserved in 
closed problems such that 

Converge /PToh C^onverge Converge Converge 

Domain Domain Domain Domain A\ 

TTTTl -1 Initial Initial Initial 

^. ^. '^'''>'^''^= E = ^total,real 

Domain Domain Domain 

which satisfies the important definite condition for closed problems that the 
total number X^Domafn^ ^^^^ molecules represented by the simulated 

molecules after convergence is equal to the total number A^totai,reai of real 
molecules in the closed physical problem {note: total here means the sum- 
mation ^Doniain ^^^^ fl^^ domain). 
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Now, we explain why the ceU's variables ntr^k'>^tr,k'>Ttr,k are updated by 
the auto-regulation schemes of Eq. (10) rather than Eq. ([3]). As we can see, 
Fi is updated smoothly by Eq. ([s]) and so the increment A^^^A^^ calculated 
by Eq. ^ is also smooth, which implies that the summations ^A/^A/, 
^(Aj^A/Q), ^(Aj^A/mcf /2) used in Eq. ( [lO| ) contain low stochastic noise. 
But, the summations ^ A/, ^(A^q), ^(A^mcf/2) over transient values as 
in Eq. ^ still have large stochastic noise due to the discontinuous events of 
simulated molecules moving into and out of cell k. 

The DSBGK algorithm described here is valid for any cell division using 
parallelepiped or tetrahedron. In the DSMC simulation of problems with 
complex configuration, we prefer to use the regular parallelepiped to divide 
the fiow domain as in p^, which makes it efficient to determine which cell 
the simulated molecules are located inside at the end of each At. Although 
the use of parallelepiped makes it time-consuming to determine the molecu- 
lar reffection position on the complex wall surface, the number of simulated 
molecules running into the surface during each At is usually much smaller 
than the total number when Kn is much smaller than 1. Compared to us- 
ing tetrahedrons to divide the ffow domain which makes the determination of 
surface reffection positions of few simulated molecules efl&cient but the deter- 
mination of the situated cells of all simulated molecules after each time step 
time-consuming, the gain of the algorithm of using parallelepiped outweighs 
its loss. But, in the DSBGK simulation, the efficiency of the algorithm of 
molecular motion and intermolecular collision processes depends less on the 
cell type because the molecular trajectories are divided into segments by 
cell's interfaces and the molecular variables are updated along each segment 
in sequence. If molecular reffections on the wall are very frequent and com- 
plex wall conffgurations are involved, we suggest to use tetrahedron to divide 
the ffow domain in the DSBGK simulation such that the determination of 
surface reffection positions is efficient. 

In the DSMC simulation, the total CPU time is almost proportional to 
the product of sample size ngampie and sampling interval ^sample as the CPU 
time used for the transitional period before reaching the steady state is usu- 
ally negligible. The molecular quantities of interest are sampled at intervals 
(<^sampie = 4 At for iustancc) to reduce the sample size risampie- We use no- 
tations VI = y (CPUtime, rfsampie,i) and V2 = y^CPUtime, rfsampie,2) to rep- 
resent the variances using different rfgampie but the same CPU time, namely 

^sample,! ^sample,! ^sample, 2 ^sample, 2- Let rfsample,2 ^ ^sample,! and SO 
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^sampie,2 < ^sample,! • Then, the general rule ^ is that 1 < ^ < ^''""'^^^^^ 

^sample, 2 

and the ratio of variance approaches to 1 when the correlation degree of sam- 
ple set is very high. So, the increase of statistical variance due to the increase 
of ^sample froui 1 to 4 uudcr the conditions of same CPU time is negligible 
because the correlation degree of consecutive transient results in the DSMC 
simulation is high. In the DSBGK simulation, the stochastic error is low 
and the sample size required to obtain smooth results is small. We prefer 
to sample ntr,^, t?tr,/c5 rtr,/c at every time step (^sample = 1) as the variance al- 
ways (although maybe slightly due to high correlation degree) decreases with 
the increase of sample size under the conditions of same CPU time due to 

1 < Yyj- Note that -yj approaches to ^^^"^p^^^^ jf the consecutive samples are 

Vl Vi ^sample,2 

almost independent, which means that the variance is inversely proportional 
to the sample size and independent of the sampling interval. 

3.3. External body force 

When considering external body force, the BGK equation is changed to: 



+ + = ^(/eq - /) (15) 



dt dxj dcj 



df , df , df 

= yueci-j) [ 

df 

where a = a{t^ x) is the acceleration due to external body force. We split — 
• , n 9f / 

mtO -^Imove ~^jQ^^ '''''^^ '^vM " / ) ^ud "^Iforce 

plify the algorithm, we decouple the effect due to — force from the other two 

ot 

effects. At the end of each At of the above DSBGK algorithm, the effects due 
df, df, 

to ^77 move 9.nd 

coll are already incorporated into the simulation and so we 

ot ot 
df, 

consider — force by changing q of each simulated molecule to q + Ata and 
_dt 

keeping x^, F/, A^^ unchanged as f{t+At^ x+ Ate, c+Ata)=f{t^ x, c) if neglect- 
ing intermolecular collision. Correspondingly, 'Utr,^ of each cell is changed to 
Utj-^k + ^ta and ntr,^, Ttr^/e keep unchanged. When sampling and outputting 
the cell's velocity, we use the average value before and after implementing a, 
namely i^tr k + 0.5Ata. 
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3.4' Boundary conditions 

For the open boundary, simulated molecules are removed from the flow 
domain when moving across the open boundary during each At. Correspond- 
ingly, some new simulated molecules are generated at the end of each At at 
the open boundary with xi and q being selected randomly as in DSMC sim- 
ulations. Then, Fi is determined from x^,q through /eq,tr,/c using the macro 
quantities fixed at the open boundary or the values of adjacent cell if not pre- 
scribed at the boundary. The initial values of Ni of new simulated molecules 
at different parts of the open boundary can be different in the DSBGK simu- 
lation. In the channel flow problem driven by the density difference Angnd at 
the two ends [3j , we use different initial values of A^^,init,end for different ends 
such that their ratios of A^^,init,end/^end are equal, which makes the number 
of simulated (not real) molecules per cell almost the same for different cells 
having the same volume but different number density of real molecules. As 
the stochastic noise at each cell depends on the average number of simu- 
lated molecules inside that cell, such selection of the initial values of A^^ for 
new simulated molecules at different parts of the open boundary achieves 
the trade-off of stochastic noise among cells and so reduces the sample size 
required for getting smooth results in the whole flow domain. 

For the wall boundary, q and then Fi are changed after molecular reflec- 
tion at xi on the wall as discussed below. 

3.4-1' Updating ci 

When running into the wall and reflecting at Xi on the wall, q is changed 
to Q,new = Cj. + tTwaii whcrc tTwaii IS the Wall velocity and the reflecting velocity 
Cj. is selected randomly according to the reflection model (namely the scatter 



kernel discussed later in section 3.4.2) as in the DSMC simulation. Ni re- 



mains unchanged to conserve the mass. After changing q alone, [x^, q, A/]aii 
is updated to represent / after molecular reflection effect and consequently 
Fi is updated to the representative value of / at the point (t, x/, Q^new)- So, 
the compatibility condition is satisfled in the molecular reflection process. 
The subscript / is omitted in the component expression of velocity when dis- 
cussing the boundary condition. We predetermine a local Cartesian reference 
system S'locai moving at the wall velocity tTwaii- We use the subscripts 2 and 
3 for the tangential directions and 1 for the normal direction of /Siocai- Iii 
the discussion of reflection process, the subscripts 1, 2, 3 always represent 
the corresponding components in /Siocai- The incoming velocity q is q — iXwaii 
{note: Q,new — '^waii = c^) • As Q and t?waii are stored in the component form 
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of the unique global Cartesian reference system S'giobab we need the transfor- 
mation from S'giobai to S'locai to obtain the components of q^i, Ci^2, Ci,3- Finally, 
Cr,i, Cr,2, Cr,3 are transformed from S'locai to ^Sgiobai to obtain the component 
form of Q^new l^ S'giobai- For the unit normal vector of wall, we specify 
that the selection of makes the incoming component q^i = q • e n negative 
and the reflecting component Cr,i = Cj. • positive. The normal direction is 
unique and the selections of tangential directions are free but flxed during 
the simulation process. In the original CLL reflection model [5j-[6j, we com- 
pute the tangential components of Cr by Cr,2 = v cos 9 — wsinO and Cr,3 = 
V8'm9 + w cos 9 where v = [{1 — o^r)(cf^2 + ^Is)]^^'^ + (2fcBrwaii/^)'^^^^r cos (fr, 
w = (2/cBTwaii/^)'^^^rT- sin^^T-, r^- = (— a^- In Lp^ = 27ri?/2, 9 is the 
azimuthal angle of incoming velocity component (ci^2,Ci,3) in the X2X3 plane 
of S\ocaii Rfi and Rf2 are two different random fractions distributed uni- 
formly inside [0, 1], ar is the accommodation coeflicient of kinetic energy 
of the tangential velocity component. For the normal component, Cr,i = 
[(1 - a^)cl + (2fcBT^aii/m)r2 + 2(1 - ,|(2fcBT^aii/m)i/V, cos (^,] 

where |ci^i| is the absolute value of q^i as q^i < 0, Tn = (— o^n In i^/a)-^/^, 
(/^n = 27ri?/4, i?/3, Rf4 are two additional random fractions and is the ac- 
commodation coeflicient of kinetic energy of the normal velocity component. 

We get Cr,2 = Ci,2(l - C^rY^'^ + (2/CBrwall/m) ^/^r^ COs((/9^ + 9) and Cr,3 = 

Ci^3(l — a^y/'^ + {2kBT^ai\/my^'^rr sin((/p^ + 9) after reorganizing the formulas 
of Cr,2, Cr,3. Notc that (f^ IS sclcctcd uniformly from a periodic interval [0, 27r] 
and so Lpr+9 can be replaced simply by (fr^ which implies that the calculation 
of 9 can be avoided to slightly improve the efliciency. So, for the CLL reflec- 
tion model, the equivalent but simpler algorithm to compute the tangential 
components in S\ocai is that = Q,2(l — ci^r)"^^^ + {"^kBT^aw/my^'^rr cos ip^ 
and Cr,3 = Ci^s{l — ary/'^ + {2kBT^ai\/my/'^rrSmipr [3j. This simpler algorithm 
also degenerates to the Maxwell diffuse reflection model when ar = = I. 

34.2. Updating Fi 

After getting Cr, Fi is updated correspondingly to F^^new = f{t) Q,new) = 
/(t, X/, Cr + tTwaii). Note that Fi is the representative value of / which is dif- 
ferent from the scatter kernel R used to select for each particular reflection 
process. Generally speaking, / is related to the mass flux but R has nothing 
to do with the mass flux. Usually, R describes the distribution probability of 
Cr inside the half velocity space {c^-e^ > 0) as a function depending on the wall 
temperature T^aii, accommodation coefficients a^^ and possibly also on the 
incoming velocity q. So, we have R = R{c^^ q) which contains T^aii, <^n, as 
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parameters. R satisfies the normalization condition ^^R{cj.^C[)dCj. = 1 
where i?(cr, Ci)dcr is the probability for the molecule coming at q to reflect 
into the velocity space element dcj. around Cj.. The transformation between 
/ at the boundary and R can be completed using the incoming mass flux. 

We introduce /b(c) as the equivalent distribution function of / observed 
in S'locai at the reflection point Xi and at the current moment t, which means 
/b(c) = /(t, xi, c + ^waii). After getting the formula of /b(c), F^^new = /B(cr). 
The distribution /B(ci)|^..^^<o of the incoming molecules is known from the 
molecular information in the adjacent cell. fB{cT)\cr-en>o is the distribution 
of reflecting molecules and related to R as introduced in 



/B(cr)(c; • e^)dc, = - R{c,,Ci)fB{ci){ci • e^)dcidc, (16) 

J Ci-en<0 



Taking integration of Eq. (16) with respect to Cj. over its half velocity space 



and using the normalization condition of i?(cr, q), we get: 

/ /B(Cr)(c; • en)dc; 

Jcr-en>0 

= - / R{c,, Ci)/B(ci)(ci • e^)dcidc, (17) 

J Cr-en>0 J Ci-en<0 

/B(Ci)(Ci • en)dci 

Ci-en<0 

which represents the mass conservation of molecular reflection process. 

777/ 

In the Maxwell diflFuse reflection model, /b difruse(cr) = neff( — — )^/^ 

-o 
— TflC 

exp( _^ ) where the eflFective n^^ will be determined by /B(ci). We as- 

2fcB-7wall 

, , .-.X / rn r-m(Ci - (2itr,/c - ^wall))^ , 

sume that /B(ci) = ntr,fe( ^ , ^ )^/^exp[ -— J where 

'^tr,k'>^tr,k'>Ttr,k ^rc the quantities of cell k close to the reflection point x/. 
Then, the number A^in of incoming real molecules on per unit wall surface 
during per unit time is: 



Nin = - /B(ci)(ci • en)dci 

J cven<0 



"'•""f^ (18) 

ntv,ky ^^^^ [exp(--»'f^) + y/7vu'ir,il + erf(M'in))] 
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where u'\^ = ^ . Similarly, the number A^out of reflecting 

y/2kBTtr,k/m 

real molecules is: 



A^out = / /B,diffuse(Cr)(c; • e^)dc, = Tl^^J (19) 



Let A^out = ^in as required in Eq. (17), we get: 



^eff = ^tr,fcA/ ^^[exp(-2i'f^) + ^/^lu'i^{l + erf(^'in))] (20) 

V ^wall 

Now, we can compute F/^new = /B,diffuse(cr) after getting n^^. We store n^^ 
and use it repeatedly for different simulated molecules reflecting on the same 
subarea close to cell k during the same At and update nefr at the end of 
each At. Additionally, the scatter kernel i?diffuse of the Maxwell diffuse re- 
flection model can be determined from /B,diffuse(cr) as we assume that it is 
independent of the incoming velocity q, namely i?diffuse = ^diffuse (cr). Using 



the formula of /B,diffuse(cr) and Eqs. ( |16| ), (17), ( |19| ) we have: 

/b ,diffuse 



-f^diffuse(Cr) 



-4e„<o/B(Ci)(Ci-e„)dCi 

/b ,diffuse 
/c..en>0 /B,diffuse(Cr)(c; • e^)&C, 

-)'exp(— — ^) 



(21) 



27r /CBTwall 2/CBrwalf 

which implies the selecting algorithm of for the Maxwell diffuse reflection 
model described in section 13.4.11 

In the CL reflection model [5j, the scatter kernel is the product of three 
independent parts related respectively to the three components: 



i?CL(cr,Ci) =^=exp[^^ ^Jx 

exp[^ X 



Cr,l r-(c?,l + (l-«n)c,?i) 

— — exp[ X 

f^^ 2 VI - ^^nCr,l|Ci,i| 

/ exp [ ■ — cos d\ ad 



(22) 
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where |ci^i| is the absolute value of the normalized incoming component 
where q i < 0. The selecting algorithm of Cj. was proposed 

y^2fcBTwaii/m 



in [6j based on Eq. (22) and is referred to as CLL reflection model. Again, 



we assume that /B(ci) is a Maxwell distribution, which is a rough assump- 



tion here although it is reasonable when calculating N^-^ by Eq. (18). Then, 



/B,CL(cr) is determined from Eqs. (16) and (22). Unfortunately, it is compli- 



cated to calculate /B,CL(cr) by solving the integral of Eq. (16). A tentative 
scheme was proposed in [2j to simplify the calculation. Note that the major 
differences between /B,difruse(cr) and i?diffuse are that the former contains a 
parameter rieff but the later contains c^-e^. So, Cr,i is removed from Rcl and 
a new parameter a is added to describe /B,CL(cr). The mass conservation 



principle of Eq. (17) is used to determine the parameter a. Consequently, 



the tentative formula of /b,cl depends not only on Cj. but also on q, which 



is inconsistent with the deflnition of Eq. (16) where /b,cl is independent 
of Q. Some simulation results show that the tentative formula of /b,cl is 
useful when a^^a^ are very close to 1 {a^ = o^n = 0.98) [3j. This is because 
the tentative formula of /b,cl can degenerate to the correct /B,diffuse when 
= = 1 and so its error is negligible when a^, o^n 1- 
In the specular reflection model, Cr = q — 2e^{ci • e^) and so i?specuiar = 

S{C, - (Ci - 2e;(Ci • en))). Submitting i?specular into Eq. (|l6|), /B,specular(Cr) = 

/B(ci), which implies F/,new = Fi as F/^new = /B,specuiar(cr) and /B(Ci) = 
f{t^xi^ci) and Fi is equal to f{t^xi^ci) before reflecting. 

3.5. Calculation of flux on boundary 

As in the DSMC method, it is convenient for the DSBGK method to 
calculate the flux T{Q) of any molecular quantity Q{c) in unit time and 
across unit area of the boundary surface: 

r(g) = (23) 

where the summation is over all those simulated molecules reflecting on the 
subarea AS during the time step At, Q{ci) and Q{cj.) are the incoming and 
reflecting quantities, respectively. Let Q = mc and mc^/2 and then T{Q) 
represents the stress and heat flux, respectively. 

3.6. Summary of the DSBGK algorithm 

1. Initialization. Generate many cells and simulated molecules and assign 
them with initial values for ntr,/c5 '^tr,/c5 ?tr,/c and x/, q, F/, A^^, respectively. 
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2. Each simulated molecule moves uniformly and in a straight line before 
encountering boundary. During each At, the trajectory of any particular 
molecule / may be divided into several segments (see Fig. [T]). Then, x^, F^, A^^ 
are updated deterministically along each segment in sequence. When encoun- 
tering the wall boundary, q is updated randomly according to the reflection 
model and then Fi is updated correspondingly. In open problems, simulated 
molecules are removed from the flow domain when moving across the open 
boundary during each At and new simulated molecules are generated at the 
open boundary at the end of each At. The variables rit^^k) ^tr, k^Tt^.k of each 
cell k is updated at the end of each At. 

3. After convergence, ntr,^, 'Str,/c5 ^tr,/c are used as the discrete solutions of 
n, 1?, T at each cell k. 



4. Simulation Results 

In DSBGK simulations, the parameter v is selected to satisfy the coef- 
ficient of viscosity /x or heat conduction k by Eq. ([6]). For problems where 
the momentum exchange is the dominant effect, we use v = v{ii) = nk^T / ji 
to satisfy ji. For problems where the heat conduction is the dominant ef- 
fect, we select v = v{k) = 5nk^T / {2mtv) to satisfy tv. Note that fi is given 
usually. For monoatomic gas where the Prandtl number is Pr = 2/3 and 
the specific heat capacity at constant pressure is Cp = 5/cB/(2m), we have 
v{hi) = 2n/cBr/(3/x) as = Cp/a/Pr. 

Lid- driven cavity flow 



y 



Argon gas: 

7o=273.15 K, r?o=2.6847xl025 ^^-3 
/uo«63 nm, Af=4.3xl0-ii s 
Tvjsw^Tq, Kn=Ao/L, adjust L=W 



Figure 2: Schematic model of lid-driven cavity fiow. 
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The results were reported first in [2]. The sizes L = W are regulated 
to change the Kn number. The Maxwell boundary condition is used and 
the cell number is 20 x 20 for Kn = 0.063 and 6.3. We set v = v{fi) in 
the DSBGK simulations. In order to reduce the infiuence due to fiuctuation 
in the number density distribution observed in the DSBGK simulations of 
closed problems (see the following Fig. [5]), the number of simulated molecules 
per cell is about 2000 in the DSBGK simulations at Kn = 0.063 and 6.3. 



DSBGK simulation at 600tli time step DSBGK simulation at 600th time step 

2000 simulated molecules per cell (no average} 2000 simulated molecules per cell (no average} 




X X 



DSBGK simulation at 600th time step DSBGK simulation at 600th time step 

2000 simulated molecules per cell (no average} 2000 simulated molecules per cell (no average} 




X X 



Figure 3: Transient results of DSBGK simulation of lid-driven problem, 
Kn = 0.063 and t/waii=0.1 m/s, 7 minutes of CPU time. 

To show the high efficiency of DSBGK simulations at low velocity, we 
choose a very small driven velocity [/waii = 0.1 m/s. Fig. [3] shows the 
transient results (no average) of DSBGK simulation at 600^^ At taking about 
7 minutes of computational time of one CPU on Lenovo laptop E43A. We 
can output many transient results at different moments of interest at the 
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additional cost of negligible computational time which is used for writing data 
to the hard disc. From the efficiency point of view, the DSBGK method is 
a promising tool for studying transient problems. But, the time coordinate 
in DSBGK simulations is not synchronous with the real time in physical 
problems due to the hysteresis effect [3j of DSBGK simulations which are 



based on the auto-regulation schemes of Eq. (10). New techniques, like time 



rescaling, are required to reduce the magnitude of hysteresis effect. 



Solid line: DSMC (time average) 
Dashed line: DSBGK (no average) 




Solid line: DSMC (time average) 
Dashed line: DSBGK (no average) 




Solid line: DSMC (time average) 
Dashed line: DSBGK (no average) 




Solid line: DSMC (time average) 
Dashed line: DSBGK (no average) 
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Figure 4: Comparison between DSMC and DSBGK methods in lid-driven 
problem, Kn 0.063 and t/waii=20 m/s. 



The driven velocity U^aW increases to 20 m/s and the transient DSBGK 
results at 600^^ At are given in Fig. [ijwith verification by the DSMC results. 
The DSBGK simulation uses about 7 minutes again, which implies that the 
computational time used by the DSBGK simulation is almost independent of 
the magnitude of deviation from equilibrium state as the average process is 
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avoided here. The DSMC simulation takes about 30 hours using 67 molecules 
per cell and about 3.4 x 10^ samples (sampling once every 4At). The compu- 
tational time required by DSMC simulation for the above case of t/waii = 0.1 
m/s can be estimated by considering the fact that the computational time 
of DSMC simulation is almost inversely proportional to the square of Mach 
number, roughly 30 x 200^ hours. 



Solid line: DSMC (time average} 

Dashed line: DSBGK at 40th time step (no average} 

Dotted line: DSBGK at 900th time step (no average} 




Solid line: DSMC (time average} 

Dashed line: DSBGK at 40th time step (no average} 

Dotted line: DSBGK at 900th time step (no average} 




Solid line: DSMC (time average} 

Dashed line: DSBGK at 40th time step (no average} 

Dotted line: DSBGK at 900th time step (no average} 




Solid line: DSMC (time average} 

Dashed line: DSBGK at 40th time step (no average} 

Dotted line: DSBGK at 900th time step (no average} 




Figure 5: Comparison between DSMC and DSBGK methods in lid-driven 
problem, Kn 6.3 and t/waii=20 m/s. 



We set [/wall = 20 m/s and increase Kn to 6.3. The DSBGK transient 
results at 40^^ At agree very well with the DSMC results. The DSBGK dis- 
tributions of T^u^v remain unchanged after 40At. But, the DSBGK dis- 
tribution of n can not stay at steady state and its deviation at 900^^ At 
from the DSMC result is remarkable. This drawback of the DSBGK method 
in closed problems implies that the ensemble- average process (if necessary) 



25 



should be used for quantities related to n. In open problems, the unphysical 
fluctuation of n is eliminated by the fixed n at open boundaries and so the 
more-efficient time-average process can be used (see the results of channel 
flow). The DSBGK simulation within 40 At takes about 11 minutes and the 
computational time for each At is increased compared to that of Kn = 0.063, 
which is because the molecular reflection on the wall becomes frequent and 
more computational time is used to generate random fractions. 

4^2. Couette flow 

The same results were reported in [3]. The distance L between two plates 
is regulated to change Kn. The cell number is 200, 20, 20 for Kn = 0.01, 0.1, 
1, respectively. The Maxwell boundary condition is used, v = v{fi) and each 
cell contains about 2550 simulated molecules in the DSBGK simulations. 



Argon gas: 

T"o=273.15 K, r)0=2.6847xl02^ m'^ 

7"waiFTo, Kn=/^/L, adjust L 
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Figure 6: Comparison between DSMC and DSBGK methods in Couette flow 
problem [3]. 
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4-3. Thermal transpiration flow 
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Figure 7: Schematic model of thermal transpiration flow. 



This problem was studied flrst in [19j where Twaii,2/rwaii,i=2. We set 
Twaii,2/rwaii,i = 1 .05ro/0.95ro ^ 1.105 to show the high efliciency of DSBGK 
simulations at low velocity. The sizes L = W are regulated to change Kn. 
The cell number is 40 x 40 for Kn = 0.2 and the Maxwell boundary condition 
is used. Each cell contains about 1000 simulated molecules and v = v{tv) in 
DSBGK simulation as the heat conduction is the dominant effect. 

The DSBGK simulation converges after 160 At taking about 8 minutes 
of computational time. The transient DSBGK results are given in Fig. |8| 
The transient n and T are smooth but the transient u and v contain large 
stochastic noise, which is because that the variation of T is the active factor 
and has strong correlation with the variation of n through the rough balance 
of pressure. However, the variations of u and v are passive factors and 
so sensitive to stochastic noise. In order to present smooth results of u 
and V for clear veriflcation, we use the time-average process after IGOAt to 
reduce noise and collect 1500 samples (sampling at each At) making the total 
computational time about 79 minutes. The DSBGK smooth results are given 
in Fig. |9]with comparison by the DSMC time-average results. The DSBGK 
results using v = v{fi) are given together to show the dependence of v on 
different problems. The comparison shows that we should select v = v{hi) in 
the thermal transpiration problem. In addition. Fig. [9] shows the agreement 
between the ensemble-average and time-average for sampling u and v in the 
DSBGK simulation, which is consistent with the conclusion of Fig. [5] that 
the time- average process is valid for sampling u and v. 
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DSBGK simulation at 160tli time step v—v(k) 
1000 simulated molecules per cell (no average) 
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DSBGK simulation at 160th time step v=v(k) 
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Figure 8: Transient results of DSBGK simulation of thermal transpiration 
problem, Kn = 0.2, r^aii,i=0.95ro, Twaii,2=1.05ro, 8 minutes of CPU time. 
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Solid line: DSMC (time average} 
Dashed line; DSBGK (no average} v=v{k) 
IE 06 Dotted line: DSBGK (no average} iy=v{pi) 
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Solid line: DSMC (time average} 

Dashed line: DSBGK (time average} v=v{k) 

Dotted line: DSBGK (time average} d=d(ju) 




4E-07 6E-07 
X 



8E-07 1 E-06 



Solid line: DSBGK (ensemble average} v=v(k) 
Dashed line: DSBGK (time average} v=d(k) 




Solid line: DSMC (time average} 
Dashed line: DSBGK (no average} o=v(k) 
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Solid line: DSBGK (ensemble average} v=v{k) 
Dashed line: DSBGK (time average} v=v{k) 




Figure 9: Comparison between DSMC and DSBGK methods in thermal 
transpiration problem, Kn = 0.2, rwaii,i=0.95To and Twaii,2=l-05To. 
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4-4' Channel flow 



Argon gas: 

T"o=273.15 K, no=2.6847xlo25 m"^ 
/.o:^63 nm, At=4.3xl0-ii s 
y T"waii=7"o. Kn=Ao/W, adjust 14/ 



/^outlet/ To 



W 



Figure 10: Schematic model of channel flow. 
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DSBGK simulation at 30000th time step 
10 simulated molecules per cell (no average) 




DSBGK simulation at 30000th time step 
10 simulated molecules per cell (no average) 




Figure 11: Transient results of DSBGK simulation of channel flow problem, 
Kn = 0.63 and noutiet=0.6no, 36 minutes of CPU time. 



The DSBGK simulations of channel flows driven by pressure difference 
were reported flrst in [3j. Here, L = 5 microns and W is regulated to change 
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Kn. The cell number is 200 x 20 for Kn = 0.63 and the Maxwell boundary 
condition is used. We set v = v{ii) in the DSBGK simulations. To show the 
stability improvement of DSBGK simulations in open problems, we appropri- 
ately choose the initial value A^/,t=o of Ni of all simulated molecules such that 
the number of simulated molecules per cell is about 10 at the initial state. 
The number density at the outlet is 0.6no and equal to the initial value no at 
the inlet. In order to make the number of simulated molecules per cell almost 
uniform during the simulation process, the initial value of Ni for the new sim- 
ulated molecules at the inlet is larger than that at the outlet and the ratio is 

A^/,init, inlet /A^/,init,outlet = ^inlet /^outlet = 1/0.6. Specifically, WC SCt A^/,init, inlet = 

A^/,f=o^iniet/^o = A^/,t=o and A^/,init,outiet = A^/,f=o^outiet/^o = 0.6A^/,t^o to main- 
tain the number of simulated molecules per cell approximately equal to 10 
during the simulation process. 

After convergence, the transient DSBGK results at 30000*'^ At are given 
in Fig. [TT] taking about 36 minutes of computational time. It shows that the 
DSBGK simulation is stable when using only 10 simulated molecules per cell 
in open problem. Using few simulated molecules reduces the memory usage 
and improves the applicability in problems of large scale. The transient n 
is smooth but the transient T, v contain large stochastic noise. We use 
the time-average process to reduce noise and collect 6000 samples (sampling 
at each At) after 30000At, which takes about 8 minutes making the total 
computational time about 44 minutes. The time-average results of DSBGK 
simulation are given in Fig. [12] with comparison by the time-average results 
of DSMC simulation. Unfortunately, the average results of v and T of the 
DSBGK and DSMC simulations are still dominated by the stochastic noise 
due to small variations inside the fiow domain, particularly in the area far 
away from the two ends. It should be pointed out that the average v and 
T can distinctly show their main variations near the inlet and outlet. Note 
that the dominance of stochastic noise is due to not only small characteristic 
velocity but also small variation. As we can see from the lid-driven problem 
at small driven velocity [/waii = 0.1 m/s, the transient velocity distribution is 
smooth during the whole evolution process as its variation inside the whole 
flow domain is obvious (see Fig. |3]). The magnitude of stochastic noise in the 
DSBGK time-average results using a small sample size is much smaller than 
that in the DSMC time-average results using a large sample size. In addi- 
tion, the agreement between the DSBGK time-average and DSBGK transient 
results of n implies that the nonphysical fluctuation of n observed in the DS- 
BGK simulation of closed problem is eliminated in the open problem and the 
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time-average process is valid for sampling n if necessary. 



Solid line: DSMC (time average} 
Dashed line: DSBGK (time average) 
Dotted line: DSBGK (no average} 
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Figure 12: Comparison between DSMC and DSBGK methods in channel 
flow problem, Kn = 0.63 and noutiet=0.6no. 
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5. Conclusions 



The DSMC algorithm is analyzed using the importance sampling scheme 
to solve the Boltzmann equation. The DSBGK algorithm is introduced by 
theoretical analysis which shows the convergence of DSBGK method to the 
BGK equation. Many numerical results in several benchmark problems are 
listed together to show the validity and high efficiency of the DSBGK method. 

Unsolved Problems in the current DSBGK algorithm include: 1) hystere- 
sis effect in transient problems; 2) nonphysical fluctuation of density distri- 
bution in dosed problems; 3) how to get a rigorous formula of /B,CL(cr) using 
Eqs. (fTel) and (I^. 
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